# konrad.smolinski@gmail.com
# date:		   23/11/2010
# last update: 29/11/2010
#
# info:
#	we assume:
#	 Y is binary: Y in {0,1}
#	 X has K-points of support: {x_1,...,x_K}
# 	 Z has R-points of support: {z_1,...,z_R}
# 	 U = (U1,U2)
#
# dependencies:
# library: gss
# sources: rcWatchOut.R
# ---------------------------------------------------------------------------------------
# normalCopula() : copula function
rcNormalCopula <- function(u=c(0.5,0.5),mu=c(0,0),sig=cbind(c(1,0),c(0,1)) ){
	x <- qnorm(u[1],mean=mu[1],sd=sqrt(sig[1,1]))
	y <- qnorm(u[2],mean=mu[2],sd=sqrt(sig[2,2]))
	rho <- sig[1,2]/sqrt(sig[1,1]*sig[2,2])
	mdnorm <- exp(-(x^2+y^2-2*rho*x*y)/(2*(1-rho^2)) )/(2*pi*sqrt(1-rho^2))
	res <- mdnorm/(dnorm(x)*dnorm(y))
return(res)
}
# ---------------------------------------------------------------------------------------
# rcIntBaseComponents() : integrated Base Components 
# dependencies: rcNormalCopula()

rcIntBaseComponents <- function(a=c(-0.1,0.5),b=c(0.1,1),d=c(0.5,0),mu=c(0,0),xval=c(-0.5,0.5),acc=35){
	# parameterization
	sigMat <- cbind(c(1,b[1]),c(b[1],b[2]+b[1]^2))	
	
	# quadratures:	
	qs <- smolyak.quad(2,acc)
	pts <- qs$pt; ww <- qs$wt
	
	# Integrand evaluated at sparse grid points
	val <- NULL
	for(i in 1:nrow(pts)) val <- c(val,rcNormalCopula(pts[i,],mu=mu,sig=sigMat) ) 
	int <- val*ww  

	# important Part:
	# support and pie slices in R^2
	ysupp <- c(FALSE,TRUE);	K <- length(x)

	# Y=0 conditional on X=x[k]
	xy0 <- NULL
	for(i in 1:K) xy0 <- cbind(xy0, (pts[,2] < pnorm(a[1]+(a[2]+qnorm(pts[,1],sd=sqrt(b[2]+b[1]^2)))*x[i],mean=mu[2],sd=1) ))

	# ----
	baseCompPtsMat <- NULL
	xysupp <- list()
	for(i in 1:K) xysupp[[i]] <- ysupp

	baseComponents <- rcBaseComponents(K)
	nOfBaseComponents <- nrow(baseComponents)

	# this can be done much shorter and more neat 
	for(i in 1:nOfBaseComponents){

		if(baseComponents[i,1]==0) baseCompPts <- xy0[,1] 
		if(baseComponents[i,1]==1) baseCompPts <- !xy0[,1] 

		for(j in 2:K){
			if(baseComponents[i,j]==0) temp <- xy0[,j] 
			if(baseComponents[i,j]==1) temp <- !xy0[,j] 
			baseCompPts <- (baseCompPts & temp)
		}
		baseCompPtsMat <- cbind(baseCompPtsMat,baseCompPts)
	}

	# Validate
	intBaseComponents <- NULL
	for(i in 1:nOfBaseComponents)	intBaseComponents <- rbind(intBaseComponents,sum(int[baseCompPtsMat[,i]]))

return(intBaseComponents)
}

sum(rcIntBaseComponents())

# --------------------------------------------------------------------------

# quadratures
accuracy <- 35

# parameterization
a0 <- -0.1; a1 <- 0.5
b1 <- 0.1 ; b2 <- 1
d1 <- 0.5 ; d2 <- 0
mu.u1 <- 0; mu.u2 <- 0
sigMat <- cbind(c(1,b1),c(b1,b2+b1^2))

# sparse grid
qs <- smolyak.quad(2,accuracy)
pts <- qs$pt; ww <- qs$wt

# Integrand evaluated at sparse grid points
val <- NULL
for(i in 1:nrow(pts)) val <- c(val,rcNormalCopula(pts[i,],mu=c(mu.u1,mu.u2),sig=sigMat) ) 
int <- val*ww  
# check: is = 1 ?
sum(int)

# important Part:
# support and pie slices in R^2
ysupp <- c(FALSE,TRUE)
x <- c(-0.5,0.5,1)
K <- length(x)

# Y=0 conditional on X=x[k]

xy0 <- NULL
for(i in 1:K) xy0 <- cbind(xy0, (pts[,2] < pnorm(a0+(a1+qnorm(pts[,1],sd=sqrt(b2+b1^2)))*x[i],mean=mu.u2,sd=1) ))

# ----
baseCompPtsMat <- NULL
xysupp <- list()
for(i in 1:K) xysupp[[i]] <- ysupp

baseComponents <- rcBaseComponents(K)
nOfBaseComponents <- nrow(baseComponents)

# this can be done much shorter and more neat 
for(i in 1:nOfBaseComponents){

	if(baseComponents[i,1]==0) baseCompPts <- xy0[,1] 
	if(baseComponents[i,1]==1) baseCompPts <- !xy0[,1] 

	for(j in 2:K){
		if(baseComponents[i,j]==0) temp <- xy0[,j] 
		if(baseComponents[i,j]==1) temp <- !xy0[,j] 
		baseCompPts <- (baseCompPts & temp)
	}
	baseCompPtsMat <- cbind(baseCompPtsMat,baseCompPts)
}

# Validate
intBaseComponents <- NULL
for(i in 1:nOfBaseComponents)	intBaseComponents <- rbind(intBaseComponents,sum(int[baseCompPtsMat[,i]]))
# does integrate to 1?	
#sum(intVec)

# Table LHS


# DONE
# --- plots ------------------------------------------------------
plt1 <- FALSE
if(plt1){
	# original variables u=(u1,u2)
	u2 <- seq(from=-10,to=10,len=100)
	u1 <- NULL
	for(i in 1:K) u1 <- cbind(u1,a0 + (a1+u2)*x[i])

	# transformed variables v=(v1,v2)
	v2 <- pnorm(u2,mean=mu.u1,sd=sqrt(b2+b1^2))
	v1 <- NULL
	for(i in 1:K) v1 <- cbind(v1,pnorm(u1[,i],mean=mu.u2,sd=1))


	par(mfrow=c(1,2))
	plot(c(-10,10),c(-10,10),ylab=expression(u[1]),xlab=expression(u[2]),col='white')
	for(i in 1:K) points(u2,u1[,i],pch=19,cex=0.1,col=i,type='l')

	abline(h=a0,col=3,lty=3)
	abline(v=-a1,col=3,lty=3)

	plot(c(0,1),c(0,1),ylab=expression(v[1]),xlab=expression(v[2]),col='white')
	for(i in 1:K) points(v2,v1[,i],pch=19,cex=0.1,col=i,type='l')

	abline(h=pnorm(a0,mean=mu.u2,sd=1),col=3,lty=3)
	abline(v=pnorm(-a1,mean=mu.u1,sd=sqrt(b2+b1^2)),col=3,lty=3)

	points(pts,pch=19,cex=0.25,col='grey')

	for(i in 1:nOfCon) points(pts[conMat[,i],],pch=19,cex=0.25,col=i)
}

plt2 <- FALSE
if(plt2){
	library(rgl)
	points3d(pts[,1],pts[,2],val)
}


#library(combinat)
#suppx <- 1:K
#combn(suppx,1)
#expand.grid(suppx,suppx)

#rcBaseComponents(K=2)